﻿/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Marcus Cuda
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

namespace dnAnalytics.LinearAlgebra.Decomposition
{
    internal class DenseHouseholder : AbstractHouseholder
    {
        private double[] mTau;

        public DenseHouseholder(Matrix matrix)
            : base(matrix)
        {
        }

        protected override void DoCompute()
        {
            mQ = new DenseMatrix(mRows, mRows);
            mTau = new double[mColumns];
            double[] q = ((DenseMatrix)mQ).Data;
            double[] r = ((DenseMatrix)mR).Data;

            SafeNativeMethods.d_qr_factor(mRows, mColumns, r, mTau, q);
        }

        protected override void DoSolve(Matrix b, Matrix result)
        {
            if (b is DenseMatrix && result is DenseMatrix)
            {
                SafeNativeMethods.d_qr_solve(mRows, mColumns, b.Columns, ((DenseMatrix)mR).Data, ((DenseMatrix)b).Data, mTau, ((DenseMatrix)result).Data);
            }
            else if (b is DenseMatrix)
            {
                DenseMatrix tmp = new DenseMatrix(result);
                SafeNativeMethods.d_qr_solve(mRows, mColumns, b.Columns, ((DenseMatrix)mR).Data, ((DenseMatrix)b).Data, mTau, (tmp).Data);
                tmp.CopyTo(result);
            }
            else if (result is DenseMatrix)
            {
                DenseMatrix tmp = new DenseMatrix(b);
                SafeNativeMethods.d_qr_solve(mRows, mColumns, b.Columns, ((DenseMatrix)mR).Data, (tmp).Data, mTau, ((DenseMatrix)result).Data);
            }
            else
            {
                DenseMatrix tmpB = new DenseMatrix(b);
                DenseMatrix tmpResult = new DenseMatrix(result);
                SafeNativeMethods.d_qr_solve(mRows, mColumns, b.Columns, ((DenseMatrix)mR).Data, tmpB.Data, mTau, tmpResult.Data);
            }
        }

        protected override void DoSolve(Vector b, Vector result)
        {
            SafeNativeMethods.d_qr_solve(mRows, mColumns, 1, ((DenseMatrix)mR).Data, ((DenseVector)b).Data, mTau, ((DenseVector)result).Data);
        }
    }
}